/************************************************************************
*************************************************************************
*************************************************************************

simulate_systematicparam_rr.do

Runs the estimators and the tests for arbitrary production function (no longer
calibrated to GNR's Chile Industry 311), with more challenging parameters
varied systematically.

*************************************************************************
************************************************************************/


args jobid sdset

local logon 1

clear
clear matrix
set more off




//Store current state of the random number generator
global curseed = c(seed)

//Load random seeds
do "do/rngseeds`sdset'.do"

//Load globals
do "constants/simglobals.do"


if $ncores >  1 {
	set seed ${seed`jobid'}
}

	
	
//Set maximum number of iterations
global prevmaxiter = c(maxiter)
set maxiter 1000


if `logon' == 1 {
	cap log close

	cap mkdir "logs"
	local tmp = string(date("`=c(current_date)'", "D M Y"),"%tdCCYYNNDD")
	cap mkdir "logs/`tmp'"
	local tmp2 = clock("`=c(current_time)'","hms")
	log using "logs/`tmp'/`tmp2'simulate_systematicparam_rr`sdset'_`jobid'"
}



//Targets for fraction of firms constrained
global consttargets 80 70 60 50 40 30 20 10 1

//This must be fixed at 10 to get all the scenarios
global subspecs 10

//DO bootstrap?
global dobs 1


//Name of folder where we'll save the result
global savename sweeps_systematicparam_widrev

/************************************************************************
*************************************************************************
Part 1: Helper Programs
************************************************************************
************************************************************************/

//Load data generator
do "do/programs/gendatsyn.do"

//Load GNR program
do "do/programs/gnr_jpe.do"

//Load ar estimation programs
do "do/programs/ar_functions.do"

//Load additional testing programs
do "do/programs/structural_idtests.do"

//Load sharetest program
do "do/programs/sharetest.do"


/************************************************************************
Program: product
Purpose: Computes binary interaction terms
************************************************************************/
cap program drop product
program product
	syntax varlist
	
	gettoken firstvar resvars : varlist	
		
	//Take the product
	foreach var of varlist `varlist' {
		cap gen `firstvar'X`var' = `firstvar'*`var'
	}
	
	//If we're not on the last word, keep going
	if !mi(word("`resvars'",1)) {
		//Continue recursion
		product `resvars'
	}
end







/************************************************************************
Program: truemoments
Purpose: Computes true moments (must tell it the true production function)

Syntax:
name	-	"ces" "quad" (anything else = cobb douglas)
pm		-	Parameter on M (theta)
pk		-	parameter on K
pl		-	parameter on L
mu		-	elasticity of substitution
************************************************************************/
cap program drop truemoments
program truemoments
	syntax name(name=name), pm(real) pk(real) pl(real) [  mu(real 0) sc(real 0) ]
	
	/************************
	Compute marginal products
	***********************/
	
	/****** CES ********/
	if "`name'" == "ces" {
		//For CES they need to pass mu
		
		if mi("`mu'") {
			error "Must pass mu (elasticity of substitution) for a CES production function"
		}
		
		if mi("`sc'") {
			local sc = 1-`pm'
		}
		
		/* There seems to be an error here
		//Marginal product of K and L
		foreach X of varlist K L {
			local x = lower("`X'")
			gen MP_`X' = `p`x''*(1-`pm') * Y^(  (1-`pm'-`mu')/(1-`pm')  ) * exp(omega + eps)^( `mu'/(1-`pm') ) / `X'^( 1 - `mu' )
		}
		*/
		
		tempvar I
		gen `I' = `pk'* K^`mu' +   `pl' * L^`mu'
		
		foreach X of varlist K L {
			local x = lower("`X'")
			
			
			
			//This is actually the elasticity times Y/X
			gen MP_`X' =  `sc' *`p`x'' * `X'^`mu' / `I'  * Y/`X'
		}
		
		//Marginal product of M is just the cobb-douglas expression
		gen MP_M = `pm' * Y/M
		
		di "This is CES..."
		
		
	}
	
	/****** Quadratic ********/
	else if "`name'" == "quad" {
		//Marginal product for K and L is just param * Y/X
		foreach X of varlist K L {
			local x = lower("`X'")
			gen MP_`X' = `p`x'' * Y/`X'
		}
		
		//Marginal product of M is slightly different
		gen MP_M = exp(omega + eps) * K^`pk' * L^`pl' * (`pm' - 2*M)
	}
	
	/****** Cobb-Douglas ********/
	
	else {
		//Marginal product is just param * Y/X
		foreach X of varlist K L M {
			local x = lower("`X'")
			gen MP_`X' = `p`x'' * Y/`X'
		}	
	}
	
	
	
	
	/************************
	Moments: Elasticities and Sdev of marginal products
	***********************/
	
	//Elasticity = MP_X * X/Y
	foreach X of varlist K L M {
		local x = lower("`X'")
		
		//This may vary by firm
		gen elas`x'_tru 	= MP_`X' * `X'/Y 
		
		//Exclude the top and bottom 1 percent
		bysort t: cumul MP_`X', gen(MPpct)
		
		bysort t: egen sdevMP`x'_tru 	= sd(MP_`X') if MPpct > .01 & MPpct < .99
		
		//Clean up
		drop MP_`X' MPpct
	}

	
end






/************************************************************************
Program: getconstraints
Purpose: Sets a global equal to the proper values of the constraint set

Arguments:
spec	-	The specification to use

************************************************************************/
cap program drop getconstraints
program getconstraints
	args spec
	
	quietly {
	
		
		preserve
		
		local count = 0
		foreach const of numlist -1(.5)32 {
			
			gendatsyn  = 2^(`const') , klfunc("${XKL3}") rho(${rhocal}) n(${nfirmnorm}) t(${nyearnorm}) pm(${pm}) sigeps(${sigeps}) sigeta(${sigeta1}) ///
								savr1(${conrho}) savr2(${conrho2}) 
								
			truemoments ces, pm(${pm}) pk(${cespk2}) pl(${cespl2}) mu(${mu2}) sc(${sc2})
			
			collapse constrained
			
			
			
			local count = `count' + 1
			local frac`count' 	= constrained in 1
			local const`count' 	= `const'
		}

		clear
		set obs `count'
		gen frac 	= .
		gen const 	= .
		forvalues i=1/`count' {
			foreach var of varlist _all {
				replace `var' = ``var'`i'' in `i'
			}
		}
		
		

		global const`spec' ""
		local conlist "${consttargets}"
		foreach target of numlist `conlist' {
			gen tmp = abs(frac-`target'/100)
			sort tmp
			local tmp1 = const in 1
			
			global const`spec' "${const`spec'} `tmp1'"
			drop tmp
		}
		
		restore
		
	}
end




/************************************************************************
Program: gnr_wrap
Purpose: A wrapper for the gnr program (for use in bootstrapping)
************************************************************************/
cap program drop gnr_wrap
program gnr_wrap
	args spec laginst
	
	cap confirm var logtmp_y
	if _rc == 0 {
		//drop y k l m
	}
	else {
		rename (y k l m) logtmp_=
	}
	
	gen share = M/Y
	//Assume we know when the true process is linear or nonlinear
	if `spec' == 6 {
		//cap gnr y k l m "GMM" "keep" pricem nonlin
		
		if !mi("`laginst'") {
			cap gnr_jpe Y L K M share, `laginst'
		}
		else {
			cap gnr_jpe Y L K M share
		}
	}
	else {
		cap gnr_jpe Y L K M share, linear `laginst'
	} 
	
	if _rc != 0 {
	
		cap gen elasm_gnr 	= .
		cap gen elask_gnr 	= .
		cap gen elasl_gnr 	= .
		cap gen ar_gnr 		= .
	
	}

	drop share

end



/************************************************************************
Program: makebssamp
Purpose: Makes a single bootstrap dataset
************************************************************************/

cap program drop makebssamp
program makebssamp

	/******
	Step 1: Prepare the original dataset for merging
	******/
	
	//Make a household identifier that is numbered consecutively
	egen idcon = group(id)
	
	//Get the max and min value of IDs
	sum idcon
	local idmax = r(max)
	local idmin = r(min)
	
	//Get the max and min values for time
	sum t
	local tmin = r(min)
	local tmax = r(max)
	
	//Number of observations
	egen tag = tag(idcon)
	count if tag
	local nobs = r(N)
	drop tag
	
	
	//Save the source
	tempfile bssource
	save `bssource'
	
	
	/******
	Step 2: Construct framework of BS sample
	******/
	
	//Draw sample
	clear
	set obs `nobs'
	gen idcon = `idmin' + int( runiform()*(`idmax'-`idmin' + 1) )
	
	//The bootstrap household identifier
	gen id_bs = _n
	
	//Add ts
	expand `=`tmax'-`tmin'+1'
	bysort id_bs: gen t = _n + `tmin'-1
	
	tab t
	
	/******
	Step 3: Bring in data from source
	******/
	merge m:1 idcon t using `bssource'
	
	//These are households in the source that were not sampled
	drop if _m == 2
	
	//These are ts in the for which a sampled household is unobserved
	drop if _m == 1
	drop _m
	
	xtset id_bs t

end








/************************************************************************
Program: bspctiles
Purpose: Bootstraps the percentiles of the rk-statistic (and the vcov of the Lambda stat)
************************************************************************/

cap program drop bspctiles
program bspctiles
	//Number of reps, stub of the estimator (gnr or ar_dem), specification,
	//and a specifier for whether we should save or not
	args nreps spec savepath
	
	/******
	Step 0: Estimate true parameters
	******/
	preserve
	sort id t
	ar_dem
	local rhoest = _b[/rho]
	foreach inst of varlist k-mXm {
		gen ENDO2_`inst' = `inst' - `rhoest'*l.`inst'
		gen INST_`inst' = l.`inst'
	}
	
	gen yrho2 = y - `rhoest'*l.y
	
	foreach endvar of varlist ENDO2_* {
		reg `endvar' INST_* l2.m l2.k l2.kXm
		estimates store `endvar'
	}
	
	ranktest_hacked (ENDO2_*) (INST_* l2.m l2.k l2.kXm), cluster(id) wald fullrank
	local nel = rowsof(r(lab))
	
	local listel
	forvalues i=1/`nel' {
		local listel `listel' lab`i'
	}
	
	
	//ivreg2 yrho2 (ENDO2_* = INST_* l2.m l2.k l2.kXm), cluster(id)
	restore
	//gen widstat_iv = e(widstat)

	/******
	Step 1: Prepare a dataset for the bootstrapped estimates
	******/
	cap postclose bsests
	
	tempfile bsests
	postfile bsests widnull `listel' using `bsests'
	
	
	/******
	Step 2: Prepare the original dataset for bootstrapping
	******/
	tempfile original
	save `original'
	
	
	
	
	/******
	Step 3: Bootstrap
	******/
	forvalues b=1/`nreps' {
		preserve
		makebssamp
		drop id idcon
		sort id_bs t
		
		//Preliminary construction
		ar_dem
		local rhoest = _b[/rho]
		foreach inst of varlist k-mXm {
			gen ENDO2_`inst' = `inst' - `rhoest'*l.`inst'
			gen INST_`inst' = l.`inst'
		}
		
		//Rank test, etc.
		ranktest_hacked (ENDO2_*) (INST_* l2.m l2.k l2.kXm), cluster(id_bs) wald fullrank
		matrix l = r(lab)
		local passlist
		forvalues i=1/`nel' {
			local passlist `passlist' ( `=el("l",`i',1)' )
		}
		
		//Nullify
		foreach endvar of varlist ENDO2_* {
			estimates restore `endvar'
			predict tmp
			replace `endvar' = `endvar' - tmp
			drop tmp
		}
		
		gen yrho2 = y - `rhoest'*l.y
		ivreg2 yrho2 (ENDO2_* = INST_* l2.m l2.k l2.kXm), cluster(id_bs)
		local widnull = e(widstat)
		

		post bsests (`widnull') `passlist'
		restore
	}
	
	postclose bsests
	
	/******
	Step 4: Compute standard errors
	******/
	
	use `bsests', clear
	
	if !mi("`savepath'") {
		save "`savepath'", replace
	}
	
	local count 0
	foreach p of numlist 10 5 1 .5 {
		local ++count
		egen tmp = pctile(widnull), p(`=100-`p'')
		sum tmp
		local p`count' = r(mean)
		drop tmp
	}
	
	cap matrix drop labv
	corr lab*, cov
	matrix labv = r(C)
	

	
	use `original', clear
	
	foreach p of numlist 1/`count' {
		gen widstatp_`p' = `p`p''
	}
	
	matrix veclabv = vec(labv)'
	svmat veclabv


end




/************************************************************************
*************************************************************************
Part 2: Scenarios
************************************************************************
************************************************************************/



/************************************************************************
Program: setparams
Purpose: Stochastically sets the parameters of the economic environment.

Syntax:
name	-	"ces" "quad" (anything else = cobb douglas)
pm		-	Parameter on M (theta)
pk		-	parameter on K
pl		-	parameter on L
mu		-	elasticity of substitution
************************************************************************/
cap program drop setparams
program setparams

	args subspec



	//Parameters of the production function
	global pm 		= .67 + 2*(runiform()-.5)*.1
	global pl 		= .28 + 2*(runiform()-.5)*.2
	global pk 		= .11 + 2*(runiform()-.5)*.06


	//global quad 	= 10



	//CES, not very cobb-douglas
	global mu2		= .8
	
	local tmp		= (runiform()-.5)
	global cespk2	= .015 	- (`tmp'-.4)*.04
	global cespl2	= .95 	+ `tmp'*1 + .1
	global sc2		= .39 	+ (runiform()-.5)*.4

	////Standard deviations
	
	
	/*********
	Default scenario
	**********/
	global sigeta1 	= .0925324
	global sigeps 	= .2542193
	global rhocal 	= .7709392
	global conrho 	= .4920887
	global conrho2 	= .3205638/.4920887 * ${conrho}

	
	/*********
	"High" scenarios
	**********/
	if `subspec' == 1 {
		global sigeta1 	= .0925324 * 2
	}
	if `subspec' == 2 {
		global sigeps 	= .2542193 * 2
	}
	if `subspec' == 3 {
		global rhocal = .7709392 + .15
	}
	if `subspec' == 4 {
		global conrho 	= .4920887 + .2
	}
	if `subspec' == 5 {
		global conrho2 	= .3205638 * 1.3
	}
	
	/*********
	"Low" scenarios
	**********/
	if `subspec' == 6 {
		global sigeta1 	= .0925324 * 1/2
	}
	if `subspec' == 7 {
		global sigeps 	= .2542193 * 1/2
	}
	if `subspec' == 8 {
		global rhocal = .7709392 - .15
	}
	if `subspec' == 9 {
		global conrho 	= .4920887 - .2
	}
	if `subspec' == 10 {
		global conrho2 	= .3205638 / 1.3
	}
	
	
	
	

	//Number of firms
	global nfirmnorm 	2613


	//Number of years
	global nyearnorm 	7

	

	/*********************
	Production Functions
	**********************/

	/*
	The DGP program can handle production functions of the form

	X(K,L)*M^pm

	OR

	X(K,L)*(pm1*M - M^2)

	where X(K,L) is arbitrary. This works because in both cases optimal M depends on
	only X(K,L), but not K and L separately.
	*/

	global XKL1 "gen X_kl = K^${pk} * L^${pl}"
	global XKL2 "gen X_kl = (  ${cespk}*K^${mu} + ${cespl}*L^${mu} ) ^ ( ${sc}/${mu} )"
	global XKL3 "gen X_kl = (  ${cespk2}*K^${mu2} + ${cespl2}*L^${mu2} ) ^ ( ${sc2}/${mu2} )"


	/*********************
	Productivity processes
	**********************/
	global nlomproc "86.30933 - 7.081727 - 33.35803*(l.omega + 7.081727) + 4.568181*(l.omega + 7.081727)^2 - .2030311*(l.omega + 7.081727)^3 "

	
end



local nspecs 1

/************************************************************************
*************************************************************************
Part 3: Sweep
************************************************************************
************************************************************************/
local nreps = ${B}

cap mkdir "dta/simulations"
cap mkdir "dta/simulations/${savename}"
cap mkdir "dta/simulations/${savename}_bs"

if ${ncores} > 1 {
	cap mkdir "dta/simulations/${savename}/`jobid'_`sdset'"
	cap mkdir "dta/simulations/${savename}_bs/`jobid'_`sdset'"
	
	local savefol "dta/simulations/${savename}/`jobid'_`sdset'"
	local bsfol "dta/simulations/${savename}_bs/`jobid'_`sdset'"
	
	di "*******************"
	di "I am job `jobid'"
	di "*******************" _n _n _n
}
else {
	local savefol "dta/simulations/${savename}"
	local bsfol "dta/simulations/${savename}_bs"
}


/*
Estimate time to complete
How many operations?
*/

local totops = `nspecs' * (${subspecs}+1) * wordcount("${consttargets}") * `nreps'
local opscom = 0
local esttime "???"

//Loop through specifications
foreach spec of numlist 1 {

	cap mkdir "`savefol'/spec`spec'"
	cap mkdir "`bsfol'/spec`spec'"
	
	
	di _n "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
	
	//Loop through assumptions about the fraction observed
	forvalues s=0/${subspecs} {
	
		di "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"
	
		//Set Parameters
		setparams `s'
		
		cap mkdir "`bsfol'/spec`spec'/`s'"
		local bshere "`bsfol'/spec`spec'/`s'"
		
	
		//Get the constraints for this specification
		getconstraints `spec'
		local conlist "${const`spec'}"
		local ind = 1
		foreach const of numlist `conlist' {
			di "-----------------------------------------------------------------------------------------"
			forvalues b=1/`nreps' {
				//Version of the constraint to display
				local constdisp = string(`const', "%9.2f")
				di "On Specification `spec'" _col(25) "Sub-spec: `s'/${subspecs}" _col(55) "Constraint = `constdisp'" _col(80) "Rep = `b'/`nreps'" _col(100) "Est. Time: `esttime' min"
			
				quietly {
					timer on 1
			
					/********
					Generate data and compute true moments
					********/
					
					gendatsyn  = 2^(`const') , klfunc("${XKL3}") rho(${rhocal}) n(${nfirmnorm}) t(${nyearnorm}) pm(${pm}) sigeps(${sigeps}) sigeta(${sigeta1}) ///
								savr1(${conrho}) savr2(${conrho2}) 
								
					truemoments ces, pm(${pm}) pk(${cespk2}) pl(${cespl2}) mu(${mu2}) sc(${sc2})
	
					order elas* sdevMP*
					
					
					

					
					/*************************************
					Step 1: Apply Methods
					*************************************/
					

					//Construct second-order terms
					order k l m, last
					product k l m
					
			
					
				
					/*********
					Share test
					*********/
					gen logshare = log( pricem*M/Y )
					sort id t
					//sharetest logshare k l m, test( l2.( c.k##c.m ) ) lagtest deg(4)
					sharetest logshare k l m, test( l2.( c.k##c.m ) ) lagtest deg(2) lagdeg(1)
					
					local stub idtest_share
					gen `stub'_pr2	= r(r2)
					gen `stub'_F	= r(F)
					gen `stub'_pval	= r(pval)
					
					
					
					sharetest logshare k l m, test( l2.( c.k##c.m ) ) mmerr deg(2)
					
					local stub idtest_mmshare
					gen `stub'_pr2	= r(r2)
					gen `stub'_F	= r(F)
					gen `stub'_pval	= r(pval)
					
					
					/********
					WIDstat test
					********/
					widstat_iv
					widstat_iv_old
					
					/********
					Get estimates
					********/
					
					foreach estimator in ar_dem gnr {
						sort id t
						`estimator'_wrap `spec' 
					}
					

					
					//Compute standard deviation of marginal products
					//foreach method in ar_std ar_dem gnr {
					foreach method in ar_dem gnr {
						foreach var of varlist K L M {
							local sm = lower("`var'")
							gen MP_`sm' = elas`sm'_`method' * Y/`var'
							
							//Exclude the top and bottom 1 percent
							bysort t: cumul MP_`sm', gen(MPpct)
							
							bysort t: egen sdevMP`sm'_`method' = sd(MP_`sm') if MPpct > .01 & MPpct < .99
							drop MP_`sm' MPpct
						}
					}
					
					
					
		
					
					sort id t
					cap drop k l
					rename logtmp_* *
					order k l m kXk kXl kXm lXl lXm mXm, last

					
										
					if ${dobs} == 1 {
						local fc : word `ind' of $consttargets
						bspctiles 100 1 "`bshere'/`fc'_`b'.dta"
					}
						
					
					
					
					/********
					Record parameters
					********/
					foreach par in pm pl pk mu2 cespk2 cespl2 sc2 sigeta1 sigeps rhocal conrho conrho2 {
						gen par_`par' = ${`par'}
					}

					
	
					
					//keep elas* sdevMP* test* widstat widpval constrained 
					keep elas*  sdevMP* constrained  par_* idtest* widstat_* widstatp_* 
					collapse _all
					
					tempfile `b'
					save ``b''
					
					
					/********
					Estimate time remaining
					********/
					timer off 1
					timer list
					local opscom 	= `opscom'+1
					local timesofar = r(t1)
					local avgtime	= `timesofar'/`opscom'
					local timerem 	= `avgtime' * (`totops'-`opscom') / 60
					local esttime 	= string(`timerem', "%9.2f" )
					
				}
				
			
			}
			
			quietly {
				use `1'
				forvalues k=2/`nreps' {
					append using ``k''
				}
				
				
				
				//if there's a negative, change it to an underscore
				gen const 	= `const'
				local fc : word `ind' of $consttargets
				gen fraccon	= `fc'
				gen subspec = `s'

				save "`savefol'/spec`spec'/subspec`s'_const`fc'.dta", replace
				
				local ind = `ind'+1
			}
		}
	}
}


//Restore default max number of iterations for GMM
set maxiter ${prevmaxiter}


if `logon' == 1 {
	cap log close
}
